Circulant SSA

L <- 81
x <- 1:L
y <- sin(pi*x/(L+1))^2 |> sqrt() 

plot(x, y)

Вспомогательные функции

library(Rssa)

Загрузка требуемого пакета: svd
Загрузка требуемого пакета: forecast
Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

Присоединяю пакет: ‘Rssa’

Следующий объект скрыт от ‘package:stats’:

    decompose

library(signal)

Присоединяю пакет: ‘signal’

Следующий объект скрыт от ‘package:Rssa’:

    roots

Следующие объекты скрыты от ‘package:stats’:

    filter, poly

library(gsignal)

Registered S3 methods overwritten by 'gsignal':
  method         from  
  plot.grpdelay  signal
  plot.specgram  signal
  print.freqs    signal
  print.freqz    signal
  print.grpdelay signal
  print.impz     signal
  print.specgram signal

Присоединяю пакет: ‘gsignal’

Следующий объект скрыт от ‘package:signal’:

    Arma, as.Arma, as.Zpg, bartlett, bilinear, blackman, boxcar, butter, buttord, cheb1ord, chebwin, cheby1, cheby2,
    chirp, conv, decimate, ellip, ellipord, fftfilt, filter, filtfilt, fir1, fir2, flattopwin, freqs, freqs_plot,
    freqz, freqz_plot, gausswin, grpdelay, hamming, hanning, ifft, impz, interp, kaiser, kaiserord, levinson, Ma,
    medfilt1, poly, remez, resample, sftrans, sgolay, sgolayfilt, specgram, triang, unwrap, Zpg, zplane

Следующие объекты скрыты от ‘package:stats’:

    filter, gaussian, poly

source("eossa_new.R")


dftmtx <- function(n) {
  y <- stats::mvfft(diag(1, n))
  y
}

diag_averaging <- function(A){
  B <- A[nrow(A):1, ] |> Re()
  lapply(split(B, -(row(B) - col(B)) ), mean) |> as.numeric()
}

shift_vector <- function(vec) {
  last_element <- tail(vec, 1)
  vec <- vec[-length(vec)]
  shifted_vec <- c(last_element, vec)
  return(shifted_vec)
}

extend <- function(x, H){
  # Вычисление коэффициентов AR модели для дифференцированного ряда
  N <- length(x)
  p <- floor(N / 3)
  dx <- diff(x)
  # A <- ar(dx, order.max = p, method = "yule-walker")$ar
  A <- aryule(dx, p)$a
  
  # Правое расширение
  y <- x
  dy <- diff(y)
  er <- signal::filter(A, 1, dy)
  dy <- signal::filter(1, A, c(er, rep(0, H)))
  y <- y[1] + c(0, cumsum(dy))
  
  # Левое расширение
  y <- rev(y)
  dy <- diff(y)
  er <- signal::filter(A,1,dy)
  dy <- signal::filter(1,A,c(er, rep(0, H)))
  y <- y[1] + c(0, cumsum(dy))
  
  # Расширенный ряд
  xe <- rev(y)
  
  # Вывод результатов
  xe 
}

CiSSA

Подаётся на вход временной ряд, длина окна (если её нет, то она равна длине ряда + 1 пополам) и информация о том, нужно ли расширить ряд. Расширять ряд стоит при стохастическом тренде (Autoregressive extension (default). It is indicated for stationary and stochastic trend time series as well). Реализовано только Autoregressive extension.

На выходе список выдаётся список list(t_series, importance).
t_series — матрица, по столбцам которой располагаются временные ряды, отвечающие за частоты (i-1)/L, где i — номер столбца, L — длина окна.
importance — вектор, отвечающий за значимость i-ого временного ряда в разлолжении. Чем больше значение, тем больший вклад внёс i-тый временной ряд.

circulant_SSA <- function(ts, L = NULL, extend_flag = FALSE){
  time_series <- ts
  # Construct trajectory matrix
  N <- length(time_series)
  if (is.null(L)){
    L <- (N + 1)%/%2
  }
  # Проверка на расширения ряда
  if (extend_flag == FALSE){
    H <- 0
    time_series <- ts
  }
  else{
    H <- L
    time_series <- extend(ts, H)
  }
  
  X <- hankel(time_series, L)
  
  # Number of symmetric frequency pairs around 1/2
  if (L %% 2) {
    nf2 <- (L + 1) / 2 - 1
  } else {
    nf2 <- L / 2 - 1
  }
  
  # Number of frequencies <= 1/2
  nft <- nf2 + abs((L %% 2) - 2)
  
  # Decomposition
  # Estimate autocovariance     OK
  autocov <- numeric(L)
  for (m in 0:(L-1)){
    autocov[[m+1]] <- sum(time_series[1:(N-m)] * time_series[(1+m):N]) / (N-m)
  }
  
  # First row of circulant matrix
  circ_first_row <- numeric(L)
  for (m in 0:(L-1)){
    circ_first_row[[m+1]] <- (L-m)/L * autocov[[m+1]] + (m)/L * autocov[[L-m]]
  }
  
  # Build circulant matrix
  S_C <- matrix(circ_first_row, nrow = 1)
  shifted_vector <- circ_first_row
  for (i in 2:(L)) {
    shifted_vector <- shift_vector(shifted_vector)
    # S_C <- rbind(S_C, as.vector(shifted_vector))
    S_C <- rbind(as.vector(shifted_vector), S_C)
  }
  
  # Eigenvectors of circulant matrix (unitary base)
  U <- dftmtx(L)/sqrt(L)
  
  # Real eigenvectors (orthonormal base)
  U[, 1] <- Re(U[, 1])
  for (k in 1:nf2) {
    u_k <- U[, k + 1]
    U[, k + 1] <- sqrt(2) * Re(u_k)
    U[, L + 2 - (k + 1)] <- sqrt(2) * Im(u_k)
  }
  if (L %% 2 != 0) {
    U[, nft] <- Re(U[, nft])
  }
  
  # Eigenvalues of circulant matrix: estimated power spectral density
  psd <- abs(diag(t(U) %*% S_C %*% U))
  
  # Principal components
  W <- t(U) %*% X
  # Reconstruction
  # Elementary reconstructed series
  R <- matrix(0, nrow = N+2*H, ncol = L)
  for (k in 1:L) {
    R[, k] <- U[ ,k] %*% t(W[k, ]) |> diag_averaging()
  }
  
  # Grouping by frequency
  # Elementary reconstructed series by frequency
  Z <- matrix(0, nrow = N+2*H, ncol = nft)
  Z[, 1] <- R[, 1]
  # Importance of component
  imp <- numeric(nft)
  lambda_sm <- sum(psd)
  imp[1] <- psd[1]/lambda_sm
  for (k in 1:nf2) {
    Z[, k + 1] <- R[, k + 1] + R[, L + 2 - (k + 1)]
    imp[k+1] <- (psd[k+1] + psd[ L + 2 - (k + 1)])/lambda_sm
  }
  if (L %% 2 != 0) {
    Z[, nft] <- R[, nft]
    imp[nft] <- psd[nft] / lambda_sm
  }
  
  list(t_series = Z[(H+1):(N+H),],
       importance = imp,
       freq = (0:(length(imp) -1))/L
       )
}

# groups - list of frequencies
grouping_cissa <- function(cissa_res, groups){
  freq <- cissa_res$freq
  t_series <- cissa_res$t_series
  
  
  residuals <- 0
  result <- setNames(as.list(rep(0, length(groups))), names(groups))
  result_freqs <- list()
  for (i in 1:length(cissa_res$freq)){
    flag <- FALSE
    for (name in names(groups)){
      if (groups[[name]][1] <= freq[i] & freq[i] <= groups[[name]][2]){
        flag <- TRUE
        result[[name]] <- result[[name]] + t_series[, i]
        result_freqs[[name]] <- c(result_freqs[[name]], freq[i])
      }
    }
    
    if (flag == FALSE){
      residuals <- residuals + t_series[, i]
    }
  }
  
  result[["residuals"]] <- residuals
  
  return(
    list(
      t_series = result,
      freqs_by_group = result_freqs
    )
  )
  result
}

generate_ts <- function(func, n=1e3, ...){
  1:n |> func(...) |> ts()
}

f_cos <- function(x, A = 1, omega = 1/4, phi = 0){
  f_exp_mod_harm_series(x, A, alpha = 0, omega = omega, phi = phi)
}

f_sin <- function(x, A = 1, omega = 1/4, phi = 3*pi/2){
  f_exp_mod_harm_series(x, A, alpha = 0, omega = omega, phi = phi)
}

f_exp <- function(x, A = 1, alpha = 1){
  A * exp(alpha * x)
}

f_exp_cos <- function(x, A = 1, alpha = 1, omega = 1/4, phi = 0){
  f_exp_mod_harm_series(x, A, alpha, omega, phi)
}

f_const <- function(x, C = 0){
  rep(C, length(x))
}

f_exp_mod_harm_series <- function(x, A = 1, alpha = 1, omega = 1/4, phi = 0){
  A*exp(alpha*x)*cos(2*pi*omega*x + phi)
}

f_linear <- function(x, a = 1, b = 0){
  a*x + b
}
mse <- function(f_true, f_reconstructed){
   mean((f_true - f_reconstructed)^2) 
}

Ошибка при Lw in N, Kw not in N

n <- 96*2+5
L <- 96
eps <- 1/97
f_sum <- function(x){
  f_const(x, C = 1) + f_cos(x, omega = 1/12) 
}


f_const |> generate_ts(n, C = 1) |>
  plot(col = "green", ylim = c(-1, 2), ylab = "f_n")
f_cos |>
  generate_ts(n, omega = 1/12) |>
  lines(col="green")

f_sum |> generate_ts(n) |> lines(lwd = 3, col='red')
f_n <- f_sum(1:n)



c <- circulant_SSA(f_n, L = 96, extend_flag = FALSE)
r <- grouping_cissa(c,
               groups = list(
                 trend = c(0, eps),
                 sesonal = c(1/12-eps, 1/12 + eps)
               )
               )$t_series

f_C <- f_const |> generate_ts(n, C = 1)
f_c <- f_cos |> generate_ts(n, omega = 1/12)
print("Ошибки при CiSSA")

[1] "Ошибки при CiSSA"

print(paste("Ошибка при вычислении C = 1: ", mse(f_C, r$trend) |> format(scientific = TRUE, digits = 2) ))

[1] "Ошибка при вычислении C = 1:  3.2e-31"

print(paste("Ошибка при вычислении cos(pi/12): ", mse(f_c, r$sesonal) |> format(scientific = TRUE, digits = 2) ))

[1] "Ошибка при вычислении cos(pi/12):  5.1e-30"

lines(1:n, r$trend, col="blue")

lines(1:n, r$sesonal, col="blue")

f_const |> generate_ts(n, C = 1) |>
  plot(col = "green", ylim = c(-1, 2), ylab = "f_n")
f_cos |>
  generate_ts(n, omega = 1/12) |>
  lines(col="green")

f_sum |> generate_ts(n) |> lines(lwd = 3, col='red')
f_n <- f_sum(1:n)

s <- ssa(f_n, L = 96)
r <- reconstruct(s, groups=list(
  trend = 1,
  sesonal = 2:3
))


print("Ошибки при SSA")

[1] "Ошибки при SSA"

print(paste("Ошибка при вычислении C = 1: ", mse(f_C, r$trend) |> format(scientific = TRUE, digits = 2)  ))

[1] "Ошибка при вычислении C = 1:  9.6e-05"

print(paste("Ошибка при вычислении cos(pi/12): ", mse(f_c, r$sesonal) |> format(scientific = TRUE, digits = 2)))

[1] "Ошибка при вычислении cos(pi/12):  9.6e-05"

lines(1:n, r$trend)

lines(1:n, r$sesonal)

Проверка разделимости непериодических компонент + автогруппировка SSA

n <- 96*2-1
L <- 96
eps <- 1/97

C <- 1
omega_cs <- 1/12
omega_sn <- 1/24
a <- 1/100
f_sum <- function(x){
  f_const(x, C = C) +
    f_cos(x, omega = omega_cs) +
    f_exp(x, a = a) +
    f_sin(x, omega = omega_sn)
}


f_C <- f_const |> generate_ts(n, C = C)
f_c <- f_cos |> generate_ts(n, omega = omega_cs)
f_s <- f_sin |> generate_ts(n, omega = omega_sn)
f_e <- f_exp |> generate_ts(n, a = a)

f_n <- f_sum(1:n)

library(xtable)

Предупреждение: пакет ‘xtable’ был собран под R версии 4.2.3

# Шаг 2: Создание примера данных
data <- data.frame(
  Метод = c("SSA", "CiSSA"),
  e_err = c(20, 20),
  c_err = c(23, 35),
  ec_err = c(20, 20),
  sin_err = c (20, 20),
  cos_err = c(1, 1)
)


# Отрисовка ряда f_n
plot(f_n, type = "l", lwd = 3, col = 'red', ylim = c(-2, 10),
     xlab = "Время", ylab = "Значения ряда", main = "Разложение временного ряда")

# Добавление отдельных компонентов (f_C, f_c, f_e)
lines(f_C, col = "blue")  # Компонент f_C

lines(f_c, col = "blue")  # Компонент f_c
lines(f_e, col = "blue")  # Компонент f_e

lines(f_s, col = "blue")

# Легенда
legend("topleft", legend = c("Весь ряд", "Компоненты"), 
       col = c("red", "blue"), lty = 1, lwd = 3)

c <- circulant_SSA(f_n, L = L, extend_flag = TRUE)
# r <- c$t_series
r <- grouping_cissa(c,
                    groups = list(
                      # trend = c(0, 1/100),
                      trend = c(0, eps),
                      sesonal_cos = c(1/12-eps, 1/12+eps),
                      sesonal_sin = c(1/24 - eps, 1/24+eps)
                    ))$t_series

data$cos_err[2] <- mse(f_s, r$sesonal_sin) |> formatC(format = "e", digits = 1)
data$sin_err[2] <- mse(f_c, r$sesonal_cos) |> formatC(format = "e", digits = 1)
data$ec_err[2] <- mse(f_C+f_e, r$trend) |> formatC(format = "e", digits = 1)


# png("C:/Users/nik1m/Desktop/уник/6 сем/курсач/Текст работы/img/trend inseparability/CiSSA.png")  # сохранение в формате PNG

plot(1:n, f_n, type = "l", lwd=3, ylim= c(-2, 10), col="red",
     xlab = "Время", ylab = "Значения ряда", main = "CiSSA разложение временного ряда")
lines(1:n, r$trend, col = "blue")

lines(1:n, r$sesonal_sin, col = "blue")
lines(1:n, r$sesonal_cos, col = "blue")

# Легенда
legend("topleft", legend = c("Весь ряд", "Компоненты"), 
       col = c("red", "blue"), lty = 1, lwd = 3)

# dev.off()  # завершение сохранения










s <- ssa(f_n, L)
# e <- eossa_new(s, nested.groups = list(1:30), clust_type = "distance")
e <- eossa(s, 1:10, k = 7)

g_sesonal <- grouping.auto(e, base = "eigen",
                   freq.bins = list(trend = c(eps),
                                    sesonal2 = c(1/24-eps, 1/24+eps),
                                    sesonal1 = c(1/12-eps, 1/12+eps)
                                    ),
                   threshold = 0.1)


r <- Rssa::reconstruct(e, groups=c(list(exp = 1,
                                C = 2
                                ),
                             g_sesonal)
                 )

plot(wcor(e, groups = 1:24), scales = list(at = c(10, 20, 30)))

Предупреждение в wcor.ossa(e, groups = 1:24) :
  Component matrices are not F-orthogonal (max F-cor is 0.93). W-cor matrix can be irrelevant

data$c_err[1] <- mse(f_C, r$C) |> formatC(format = "e", digits = 1)
data$e_err[1] <- mse(f_e, r$exp) |> formatC(format = "e", digits = 1)
data$cos_err[1] <- mse(f_c, r$sesonal1) |> formatC(format = "e", digits = 1)
data$sin_err[1] <- mse(f_s, r$sesonal2) |> formatC(format = "e", digits = 1)
data$ec_err[1] <- mse(f_C+f_e, r$C+r$exp) |> formatC(format = "e", digits = 1)


# png("C:/Users/nik1m/Desktop/уник/6 сем/курсач/Текст работы/img/trend inseparability/SSA.png")  # сохранение в формате PNG

plot(1:n, f_n, type = "l", lwd=3, ylim= c(-2, 10), col="red",
     xlab = "Время", ylab = "Значения ряда", main = "SSA разложение временного ряда")

lines(1:n, r$trend, type = "l", col="green")

lines(1:n, r$exp, type = "l", ylim= c(-2, 10), col="blue")
lines(1:n, r$C, col = "blue")

lines(1:n, r$sesonal1, col = "blue")
lines(1:n, r$sesonal2, col = "blue")

# Легенда
legend("topleft", legend = c("Весь ряд", "Компоненты"), 
       col = c("red", "blue"), lty = 1, lwd = 3)

# Шаг 3: Преобразование данных в формат LaTeX
table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Thu Nov 28 18:25:06 2024
\begin{table}[ht]
\centering
\begin{tabular}{llllll}
  \hline
Метод & e\_err & c\_err & ec\_err & sin\_err & cos\_err \\ 
  \hline
SSA & 2.2e-25 & 2.2e-25 & 4.2e-28 & 3.8e-29 & 1.6e-29 \\ 
  CiSSA & 20 & 35 & 3.5e-02 & 7.7e-04 & 1.9e-03 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Пример cos*exp

n <- 96*2-1
L <- 96
eps <- 1/(L+1)

C <- 1
omega_cs <- 1/12
omega_sn <- 1/24
a <- 1/100
omega_exp <- 1/48
f_sum <- function(x){
    f_cos(x, omega = omega_cs) +
    f_exp_mod_harm_series(x, a = a, omega = omega_exp) +
    f_sin(x, omega = omega_sn) 
}


f_c <- f_cos |> generate_ts(n, omega = omega_cs)
f_s <- f_sin |> generate_ts(n, omega = omega_sn)
f_e <- f_exp_mod_harm_series |> generate_ts(n, a = a, omega = omega_exp)

f_n <- f_sum(1:n)

library(xtable)

# Шаг 2: Создание примера данных
data <- data.frame(
  Метод = c("SSA", "CiSSA"),
  exp_err = c(20, 20),
  sin_err = c (20, 20),
  cos_err = c(1, 1)
)


# Отрисовка ряда f_n
plot(f_n, type = "l", lwd = 3, col = 'red', ylim = c(-10, 10),
     xlab = "Время", ylab = "Значения ряда", main = "Разложение временного ряда")

# Добавление отдельных компонентов (f_C, f_c, f_e)
lines(f_c, col = "blue")  # Компонент f_c

lines(f_e, col = "blue")  # Компонент f_e
lines(f_s, col = "blue")

# Легенда
legend("topleft", legend = c("Весь ряд", "Компоненты"), 
       col = c("red", "blue"), lty = 1, lwd = 3)

c <- circulant_SSA(f_n, L = L, extend_flag = TRUE)
# r <- c$t_series
r <- grouping_cissa(c,
                    groups = list(
                      trend = c(0, 1/26-eps),
                      sesonal_cos = c(1/12-eps, 1/12+eps),
                      sesonal_sin = c(1/24-eps, 1/24+eps)
                    ))$t_series

data$cos_err[2] <- mse(f_s, r$sesonal_sin) |> formatC(format = "e", digits = 1)
data$sin_err[2] <- mse(f_c, r$sesonal_cos) |> formatC(format = "e", digits = 1)
data$exp_err[2] <- mse(f_e, r$trend) |> formatC(format = "e", digits = 1)


# png("C:/Users/nik1m/Desktop/уник/6 сем/курсач/Текст работы/img/trend inseparability/CiSSA.png")  # сохранение в формате PNG

plot(1:n, f_n, type = "l", lwd=3, ylim= c(-10, 10), col="red",
     xlab = "Время", ylab = "Значения ряда", main = "CiSSA разложение временного ряда")
lines(1:n, r$trend, col = "blue")

lines(1:n, r$sesonal_sin, col = "blue")
lines(1:n, r$sesonal_cos, col = "blue")

# Легенда
legend("topleft", legend = c("Весь ряд", "Компоненты"), 
       col = c("red", "blue"), lty = 1, lwd = 3)

# dev.off()  # завершение сохранения










s <- ssa(f_n, L)
e <- eossa_new(s, nested.groups = list(1:30), clust_type = "distance")

g_sesonal <- grouping.auto(e, base = "eigen",
                   freq.bins = list(trend = c(1/25-eps),
                                    sesonal2 = c(1/24-eps, 1/24+eps),
                                    sesonal1 = c(1/12-eps, 1/12+eps)
                                    ),
                   threshold = 0.1)


r <- reconstruct(e, groups= g_sesonal)

plot(wcor(e, groups = 1:24), scales = list(at = c(10, 20, 30)))

Предупреждение в wcor.ossa(e, groups = 1:24) :
  Component matrices are not F-orthogonal (max F-cor is -0.529). W-cor matrix can be irrelevant

data$exp_err[1] <- mse(f_e, r$trend)  |> formatC(format = "e", digits = 1)
data$cos_err[1] <- mse(f_c, r$sesonal1) |> formatC(format = "e", digits = 1)
data$sin_err[1] <- mse(f_s, r$sesonal2) |> formatC(format = "e", digits = 1)


# png("C:/Users/nik1m/Desktop/уник/6 сем/курсач/Текст работы/img/trend inseparability/SSA.png")  # сохранение в формате PNG

plot(1:n, f_n, type = "l", lwd=3, ylim= c(-10, 10), col="red",
     xlab = "Время", ylab = "Значения ряда", main = "SSA разложение временного ряда")

lines(1:n, r$trend, type = "l", col="blue")

lines(1:n, r$sesonal1, col = "blue")
lines(1:n, r$sesonal2, col = "blue")

# Легенда
legend("topleft", legend = c("Весь ряд", "Компоненты"), 
       col = c("red", "blue"), lty = 1, lwd = 3)

# Шаг 3: Преобразование данных в формат LaTeX
table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Thu Nov 28 18:25:07 2024
\begin{table}[ht]
\centering
\begin{tabular}{llll}
  \hline
Метод & exp\_err & sin\_err & cos\_err \\ 
  \hline
SSA & 4.7e-29 & 1.1e-29 & 8.4e-30 \\ 
  CiSSA & 3.2e-02 & 1.0e-03 & 5.8e-03 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Данные IP

library(readxl)

Предупреждение: пакет ‘readxl’ был собран под R версии 4.2.3

data <- read_excel("Data/International_Financial_Statistics_.xlsx")

New names:
• `` -> `...2`
• `` -> `...3`
• `` -> `...4`
• `` -> `...5`
• `` -> `...6`
• `` -> `...7`
• `` -> `...8`
• `` -> `...9`
• `` -> `...10`
• `` -> `...11`
• `` -> `...12`
• `` -> `...13`
• `` -> `...14`
• `` -> `...15`
• `` -> `...16`
• `` -> `...17`
• `` -> `...18`
• `` -> `...19`
• `` -> `...20`
• `` -> `...21`
• `` -> `...22`
• `` -> `...23`
• `` -> `...24`
• `` -> `...25`
• `` -> `...26`
• `` -> `...27`
• `` -> `...28`
• `` -> `...29`
• `` -> `...30`
• `` -> `...31`
• `` -> `...32`
• `` -> `...33`
• `` -> `...34`
• `` -> `...35`
• `` -> `...36`
• `` -> `...37`
• `` -> `...38`
• `` -> `...39`
• `` -> `...40`
• `` -> `...41`
• `` -> `...42`
• `` -> `...43`
• `` -> `...44`
• `` -> `...45`
• `` -> `...46`
• `` -> `...47`
• `` -> `...48`
• `` -> `...49`
• `` -> `...50`
• `` -> `...51`
• `` -> `...52`
• `` -> `...53`
• `` -> `...54`
• `` -> `...55`
• `` -> `...56`
• `` -> `...57`
• `` -> `...58`
• `` -> `...59`
• `` -> `...60`
• `` -> `...61`
• `` -> `...62`
• `` -> `...63`
• `` -> `...64`
• `` -> `...65`
• `` -> `...66`
• `` -> `...67`
• `` -> `...68`
• `` -> `...69`
• `` -> `...70`
• `` -> `...71`
• `` -> `...72`
• `` -> `...73`
• `` -> `...74`
• `` -> `...75`
• `` -> `...76`
• `` -> `...77`
• `` -> `...78`
• `` -> `...79`
• `` -> `...80`
• `` -> `...81`
• `` -> `...82`
• `` -> `...83`
• `` -> `...84`
• `` -> `...85`
• `` -> `...86`
• `` -> `...87`
• `` -> `...88`
• `` -> `...89`
• `` -> `...90`
• `` -> `...91`
• `` -> `...92`
• `` -> `...93`
• `` -> `...94`
• `` -> `...95`
• `` -> `...96`
• `` -> `...97`
• `` -> `...98`
• `` -> `...99`
• `` -> `...100`
• `` -> `...101`
• `` -> `...102`
• `` -> `...103`
• `` -> `...104`
• `` -> `...105`
• `` -> `...106`
• `` -> `...107`
• `` -> `...108`
• `` -> `...109`
• `` -> `...110`
• `` -> `...111`
• `` -> `...112`
• `` -> `...113`
• `` -> `...114`
• `` -> `...115`
• `` -> `...116`
• `` -> `...117`
• `` -> `...118`
• `` -> `...119`
• `` -> `...120`
• `` -> `...121`
• `` -> `...122`
• `` -> `...123`
• `` -> `...124`
• `` -> `...125`
• `` -> `...126`
• `` -> `...127`
• `` -> `...128`
• `` -> `...129`
• `` -> `...130`
• `` -> `...131`
• `` -> `...132`
• `` -> `...133`
• `` -> `...134`
• `` -> `...135`
• `` -> `...136`
• `` -> `...137`
• `` -> `...138`
• `` -> `...139`
• `` -> `...140`
• `` -> `...141`
• `` -> `...142`
• `` -> `...143`
• `` -> `...144`
• `` -> `...145`
• `` -> `...146`
• `` -> `...147`
• `` -> `...148`
• `` -> `...149`
• `` -> `...150`
• `` -> `...151`
• `` -> `...152`
• `` -> `...153`
• `` -> `...154`
• `` -> `...155`
• `` -> `...156`
• `` -> `...157`
• `` -> `...158`
• `` -> `...159`
• `` -> `...160`
• `` -> `...161`
• `` -> `...162`
• `` -> `...163`
• `` -> `...164`
• `` -> `...165`
• `` -> `...166`
• `` -> `...167`
• `` -> `...168`
• `` -> `...169`
• `` -> `...170`
• `` -> `...171`
• `` -> `...172`
• `` -> `...173`
• `` -> `...174`
• `` -> `...175`
• `` -> `...176`
• `` -> `...177`
• `` -> `...178`
• `` -> `...179`
• `` -> `...180`
• `` -> `...181`
• `` -> `...182`
• `` -> `...183`
• `` -> `...184`
• `` -> `...185`
• `` -> `...186`
• `` -> `...187`
• `` -> `...188`
• `` -> `...189`
• `` -> `...190`
• `` -> `...191`
• `` -> `...192`
• `` -> `...193`
• `` -> `...194`
• `` -> `...195`
• `` -> `...196`
• `` -> `...197`
• `` -> `...198`
• `` -> `...199`
• `` -> `...200`
• `` -> `...201`
• `` -> `...202`
• `` -> `...203`
• `` -> `...204`
• `` -> `...205`
• `` -> `...206`
• `` -> `...207`
• `` -> `...208`
• `` -> `...209`
• `` -> `...210`
• `` -> `...211`
• `` -> `...212`
• `` -> `...213`
• `` -> `...214`
• `` -> `...215`
• `` -> `...216`
• `` -> `...217`
• `` -> `...218`
• `` -> `...219`
• `` -> `...220`
• `` -> `...221`
• `` -> `...222`
• `` -> `...223`
• `` -> `...224`
• `` -> `...225`
• `` -> `...226`
• `` -> `...227`
• `` -> `...228`
• `` -> `...229`
• `` -> `...230`
• `` -> `...231`
• `` -> `...232`
• `` -> `...233`
• `` -> `...234`
• `` -> `...235`
• `` -> `...236`
• `` -> `...237`
• `` -> `...238`
• `` -> `...239`
• `` -> `...240`
• `` -> `...241`
• `` -> `...242`
• `` -> `...243`
• `` -> `...244`
• `` -> `...245`
• `` -> `...246`
• `` -> `...247`
• `` -> `...248`
• `` -> `...249`
• `` -> `...250`
• `` -> `...251`
• `` -> `...252`
• `` -> `...253`
• `` -> `...254`
• `` -> `...255`
• `` -> `...256`
• `` -> `...257`
• `` -> `...258`
• `` -> `...259`
• `` -> `...260`
• `` -> `...261`
• `` -> `...262`
• `` -> `...263`
• `` -> `...264`
• `` -> `...265`
• `` -> `...266`
• `` -> `...267`
• `` -> `...268`
• `` -> `...269`
• `` -> `...270`
• `` -> `...271`
• `` -> `...272`
• `` -> `...273`
• `` -> `...274`
• `` -> `...275`
• `` -> `...276`
• `` -> `...277`
• `` -> `...278`
• `` -> `...279`
• `` -> `...280`
• `` -> `...281`
• `` -> `...282`
• `` -> `...283`
• `` -> `...284`
• `` -> `...285`
• `` -> `...286`
• `` -> `...287`
• `` -> `...288`
• `` -> `...289`
• `` -> `...290`
• `` -> `...291`
• `` -> `...292`
• `` -> `...293`
• `` -> `...294`
• `` -> `...295`
• `` -> `...296`
• `` -> `...297`
• `` -> `...298`
• `` -> `...299`
• `` -> `...300`
• `` -> `...301`
• `` -> `...302`
• `` -> `...303`
• `` -> `...304`
• `` -> `...305`
• `` -> `...306`
• `` -> `...307`
• `` -> `...308`
• `` -> `...309`
• `` -> `...310`
• `` -> `...311`
• `` -> `...312`
• `` -> `...313`
• `` -> `...314`
• `` -> `...315`
• `` -> `...316`
• `` -> `...317`
• `` -> `...318`
• `` -> `...319`
• `` -> `...320`
• `` -> `...321`
• `` -> `...322`
• `` -> `...323`
• `` -> `...324`
• `` -> `...325`
• `` -> `...326`
• `` -> `...327`
• `` -> `...328`
• `` -> `...329`
• `` -> `...330`
• `` -> `...331`
• `` -> `...332`
• `` -> `...333`
• `` -> `...334`
• `` -> `...335`
• `` -> `...336`
• `` -> `...337`
• `` -> `...338`
• `` -> `...339`
• `` -> `...340`
• `` -> `...341`
• `` -> `...342`
• `` -> `...343`
• `` -> `...344`
• `` -> `...345`
• `` -> `...346`
• `` -> `...347`
• `` -> `...348`
• `` -> `...349`
• `` -> `...350`
• `` -> `...351`
• `` -> `...352`
• `` -> `...353`
• `` -> `...354`
• `` -> `...355`
• `` -> `...356`
• `` -> `...357`
• `` -> `...358`
• `` -> `...359`
• `` -> `...360`
• `` -> `...361`
• `` -> `...362`
• `` -> `...363`
• `` -> `...364`
• `` -> `...365`
• `` -> `...366`
• `` -> `...367`
• `` -> `...368`
• `` -> `...369`
• `` -> `...370`
• `` -> `...371`
• `` -> `...372`
• `` -> `...373`
• `` -> `...374`
• `` -> `...375`
• `` -> `...376`
• `` -> `...377`
• `` -> `...378`
• `` -> `...379`
• `` -> `...380`
• `` -> `...381`
• `` -> `...382`
• `` -> `...383`
• `` -> `...384`
• `` -> `...385`
• `` -> `...386`
• `` -> `...387`
• `` -> `...388`
• `` -> `...389`
• `` -> `...390`
• `` -> `...391`
• `` -> `...392`
• `` -> `...393`
• `` -> `...394`
• `` -> `...395`
• `` -> `...396`
• `` -> `...397`
• `` -> `...398`
• `` -> `...399`
• `` -> `...400`
• `` -> `...401`
• `` -> `...402`
• `` -> `...403`
• `` -> `...404`
• `` -> `...405`
• `` -> `...406`
• `` -> `...407`
• `` -> `...408`
• `` -> `...409`
• `` -> `...410`
• `` -> `...411`
• `` -> `...412`
• `` -> `...413`
• `` -> `...414`
• `` -> `...415`
• `` -> `...416`
• `` -> `...417`
• `` -> `...418`
• `` -> `...419`
• `` -> `...420`
• `` -> `...421`
• `` -> `...422`
• `` -> `...423`
• `` -> `...424`
• `` -> `...425`
• `` -> `...426`
• `` -> `...427`
• `` -> `...428`
• `` -> `...429`
• `` -> `...430`
• `` -> `...431`
• `` -> `...432`
• `` -> `...433`
• `` -> `...434`
• `` -> `...435`
• `` -> `...436`
• `` -> `...437`
• `` -> `...438`
• `` -> `...439`
• `` -> `...440`
• `` -> `...441`
• `` -> `...442`
• `` -> `...443`
• `` -> `...444`
• `` -> `...445`
• `` -> `...446`
• `` -> `...447`
• `` -> `...448`
• `` -> `...449`
• `` -> `...450`
• `` -> `...451`
• `` -> `...452`
• `` -> `...453`
• `` -> `...454`
• `` -> `...455`
• `` -> `...456`
• `` -> `...457`
• `` -> `...458`
• `` -> `...459`
• `` -> `...460`
• `` -> `...461`
• `` -> `...462`
• `` -> `...463`
• `` -> `...464`
• `` -> `...465`
• `` -> `...466`
• `` -> `...467`
• `` -> `...468`
• `` -> `...469`
• `` -> `...470`
• `` -> `...471`
• `` -> `...472`
• `` -> `...473`
• `` -> `...474`
• `` -> `...475`
• `` -> `...476`
• `` -> `...477`
• `` -> `...478`
• `` -> `...479`
• `` -> `...480`
• `` -> `...481`
• `` -> `...482`
• `` -> `...483`
• `` -> `...484`
• `` -> `...485`
• `` -> `...486`
• `` -> `...487`
• `` -> `...488`
• `` -> `...489`
• `` -> `...490`
• `` -> `...491`
• `` -> `...492`
• `` -> `...493`
• `` -> `...494`
• `` -> `...495`
• `` -> `...496`
• `` -> `...497`
• `` -> `...498`
• `` -> `...499`
• `` -> `...500`
• `` -> `...501`
• `` -> `...502`
• `` -> `...503`
• `` -> `...504`
• `` -> `...505`
• `` -> `...506`
• `` -> `...507`
• `` -> `...508`
• `` -> `...509`
• `` -> `...510`
• `` -> `...511`
• `` -> `...512`
• `` -> `...513`
• `` -> `...514`
• `` -> `...515`
• `` -> `...516`
• `` -> `...517`
• `` -> `...518`
• `` -> `...519`
• `` -> `...520`
• `` -> `...521`
• `` -> `...522`
• `` -> `...523`
• `` -> `...524`
• `` -> `...525`
• `` -> `...526`
• `` -> `...527`
• `` -> `...528`
• `` -> `...529`
• `` -> `...530`
• `` -> `...531`
• `` -> `...532`
• `` -> `...533`
• `` -> `...534`
• `` -> `...535`
• `` -> `...536`
• `` -> `...537`
• `` -> `...538`
• `` -> `...539`
• `` -> `...540`
• `` -> `...541`
• `` -> `...542`
• `` -> `...543`
• `` -> `...544`
• `` -> `...545`
• `` -> `...546`
• `` -> `...547`
• `` -> `...548`
• `` -> `...549`
• `` -> `...550`
• `` -> `...551`
• `` -> `...552`
• `` -> `...553`
• `` -> `...554`
• `` -> `...555`
• `` -> `...556`
• `` -> `...557`
• `` -> `...558`
• `` -> `...559`
• `` -> `...560`
• `` -> `...561`
• `` -> `...562`
• `` -> `...563`
• `` -> `...564`
• `` -> `...565`
• `` -> `...566`
• `` -> `...567`
• `` -> `...568`
• `` -> `...569`
• `` -> `...570`
• `` -> `...571`
• `` -> `...572`
• `` -> `...573`
• `` -> `...574`
• `` -> `...575`
• `` -> `...576`
• `` -> `...577`
• `` -> `...578`
• `` -> `...579`

data |> head()

Отрисовка данных IP

dates <- seq(as.Date("1970-01-01"), as.Date("2018-1-30"), by = "month")
IP_values <- data[2, -c(1, 2)] |> as.double() 
plot(dates, IP_values, type="l")

Cissa

Отрисовка трендовой составляющей чёрным цветом, основной временной ряд — красным

data_slice <- 1:537
dates_slice <- dates[data_slice]
IP_values_slice <- IP_values[data_slice]
eps <- 1/193

c <- circulant_SSA(IP_values_slice, L = 192, extend_flag = TRUE)
r <- c$t_series
r <- grouping_cissa(c,
                    groups = list(
                      trend = c(0, 1/192),
                      cycle = c(1/97, 5/95),
                      sesonal = c(1/13, 1/2+0.0001)
                    )
                    )$t_series
r_sesonal <-  grouping_cissa(c,
                             groups = list(
                              s1 = c(16/192 - eps, 16/192 + eps),
                              s2 = c(32/192 - eps, 32/192 + eps),
                              s3 = c(48/192 - eps, 48/192 + eps),
                              s4 = c(64/192 - eps, 64/192 + eps),
                              s5 = c(80/192 - eps, 80/192 + eps),
                              s6 = c(96/192 - eps, 96/192 + eps)
                             )
                             )$t_series
# cissa_trend <- r[,1] + r[,2]
# cissa_cycle <- r[, 3:11] |> rowSums()
# cissa_sesonal <- r[, c(17, 33, 49, 65, 81, 97)] |> rowSums()
# cissa_residuals <- IP_values_slice - (cissa_trend + cissa_cycle + cissa_sesonal)

cissa_trend <- r$trend
cissa_cycle <- r$cycle
cissa_sesonal <- Reduce("+", r_sesonal |> within(rm(residuals)))
cissa_residuals <- IP_values_slice - (cissa_trend + cissa_cycle + cissa_sesonal)


plot(dates_slice, IP_values_slice,
     type="l", col = "black")
lines(dates_slice, cissa_trend,
      type="l", col = "red")

plot(dates_slice, cissa_cycle,
     type="l", col = "red")

plot(dates_slice, cissa_sesonal,
     type="l", col = "red")

plot(dates_slice, cissa_residuals,
     type="l", col = "red")

plot(dates_slice, IP_values_slice,
     type="l", col = "black")
lines(dates_slice, cissa_trend+cissa_cycle+cissa_sesonal,
      type="l", col = "red")

SSA fossa

s <- ssa(IP_values_slice, L = 192)
e <- fossa(s)
# e <- eossa_new(s, nested.groups = list(1:30), clust_type = "distance")
eps <- 1/193

groups <- grouping.auto(e,
                   freq.bins = list(trend = c(1/192),
                                    cycle = c(1/97, 5/95),
                                    s1 = c(16/192 - eps, 16/192 + eps),
                                    s2 = c(32/192 - eps, 32/192 + eps),
                                    s3 = c(48/192 - eps, 48/192 + eps),
                                    s4 = c(64/192 - eps, 64/192 + eps),
                                    s5 = c(80/192 - eps, 80/192 + eps),
                                    s6 = c(96/192 - eps, 96/192 + eps)
                                    ),
                   threshold = 0)


plot(wcor(e, groups = 1:30), scales = list(at = c(10, 20, 30)),
     main = "W-correlation matrix SSA (fossa)")

r <- reconstruct(e, groups=groups)

ssa_trend_f <- r$trend
ssa_cycle_f <- r$cycle
ssa_sesonal_f <- r$s1 + r$s2 + r$s3 + r$s4 + r$s5 + r$s6
ssa_residuals_f <- IP_values_slice - (ssa_trend_f + ssa_cycle_f + ssa_sesonal_f)

plot(dates_slice, IP_values_slice,
     type="l", col = "black")
lines(dates_slice, ssa_trend_f,
      type="l", col = "magenta")

plot(dates_slice, ssa_cycle_f, 
     type="l", col = "magenta")

plot(dates_slice, ssa_sesonal_f, 
     type="l", col = "magenta")

plot(dates_slice, ssa_residuals_f,
     type="l", col = "magenta")

SSA eossa

library(Rssa)
source("eossa_new.r")
s <- ssa(IP_values_slice, L = 192)
e <- eossa_new(s, nested.groups = list(1:30), clust_type = "distance")




groups <- grouping.auto(e,
                   freq.bins = list(trend = c(1/192),
                                    cycle = c(1/97, 5/95),
                                    s1 = c(16/192 - eps, 16/192 + eps),
                                    s2 = c(32/192 - eps, 32/192 + eps),
                                    s3 = c(48/192 - eps, 48/192 + eps),
                                    s4 = c(64/192 - eps, 64/192 + eps),
                                    s5 = c(80/192 - eps, 80/192 + eps),
                                    s6 = c(96/192 - eps, 96/192 + eps)
                                    ),
                   threshold = 0)
plot(wcor(e, groups = 1:30), scales = list(at = c(10, 20, 30)),
     main = "W-correlation matrix SSA (eossa)")

Предупреждение в wcor.ossa(e, groups = 1:30) :
  Component matrices are not F-orthogonal (max F-cor is -0.0621). W-cor matrix can be irrelevant

r <- reconstruct(e, groups=groups)

ssa_trend <- r$trend
ssa_cycle <- r$cycle
ssa_sesonal <- r$s1 + r$s2 + r$s3 + r$s4 + r$s5 + r$s6
ssa_residuals <- IP_values_slice - (ssa_trend + ssa_cycle + ssa_sesonal)

plot(dates_slice, IP_values_slice,
     type="l", col = "black")
lines(dates_slice, ssa_trend,
      type="l", col = "blue")

plot(dates_slice, ssa_cycle, 
     type="l", col = "blue")

plot(dates_slice, ssa_sesonal, 
     type="l", col = "blue")

plot(dates_slice, ssa_residuals,
     type="l", col = "blue")

plot(dates_slice, IP_values_slice,
     main = "IP USA тренд",xlab = "Время", ylab = "Значение",
     type="l", col = "black")
lines(dates_slice, ssa_trend,
      type="l", col = "blue", lwd=2)

lines(dates_slice, ssa_trend_f,
      type="l", col = "magenta", lwd=2)
lines(dates_slice, cissa_trend,
      type="l", col = "red", lwd=2)

# Легенда
legend("topleft", legend = c("Весь ряд", "CiSSA тренд", "SSA тренд (eossa)", "SSA тренд (fossa)"), 
       col = c("black", "red", "blue", "magenta"), lty = 1, lwd = 3)

plot(dates_slice, ssa_cycle,
     main = "IP USA цикличность", xlab = "Время", ylab = "Значение",
     type="l", col = "blue", ylim=c(-10, 10), lwd=2)
lines(dates_slice, cissa_cycle,
      type="l", col = "red", lwd=2)

lines(dates_slice, ssa_cycle_f,
      type="l", col = "magenta", lwd=2)
# Легенда
legend("topleft", legend = c("CiSSA", "SSA (eossa)", "SSA (fossa)"), 
       col = c("red", "blue", "magenta"), lty = 1, lwd = 3)

# Настройка графиков для отображения двух графиков один под другим с общей осью X
layout(matrix(c(1, 2), nrow = 2, byrow = TRUE), heights = c(1, 1.2))

# Построение первого графика
par(mar = c(2, 4, 2, 2)) # Уменьшение нижнего отступа
plot(dates_slice, ssa_sesonal, type = "l", col = "blue", lwd = 1,
     main = "SSA (eossa) сезонность", xlab = "", ylab = "Значение")
# Добавление оси X внизу первого графика, но с пустыми метками
axis(1, labels = FALSE)

# Построение второго графика
par(mar = c(5, 4, 2, 2)) # Увеличение нижнего отступа
plot(dates_slice, ssa_sesonal_f, type = "l", col = "magenta", lwd = 1,
     main = "SSA (fossa) сезонность", xlab = "Время", ylab = "Значение")

par(mar = c(3, 4, 2, 2)) # Увеличение нижнего отступа

plot(dates_slice, cissa_sesonal, type = "l", col = "red", lwd = 1,
     main = "CiSSA сезонность", xlab = "Время", ylab = "Значение")

# Восстановление макета по умолчанию
layout(1)

NA
NA

plot(dates_slice, ssa_residuals, 
     main = "IP USA остаток", xlab = "Время", ylab = "Значение",
     type="l", col = "blue", ylim=c(-2, 2))
lines(dates_slice, cissa_residuals,
      type="l", col = "red")

lines(dates_slice, ssa_residuals_f,
      type="l", col = "magenta")
legend("topleft", legend = c("CiSSA", "SSA (eossa)", "SSA (fossa)"), 
       col = c("red", "blue", "magenta"), lty = 1, lwd = 3)

ssa_residuals |> density() |> plot()

cissa_residuals |> density() |> plot()

Отделение сигнала от шума

set.seed(100)

n_mse_tests <- function(n){
  n <- 96*2-1
  L <- 96
  sigma <- 0.1
  
  
  C <- 1
  omega_cs <- 1/12
  omega_sn <- 1/24
  a <- 1/100
  f_sum <- function(x){
    f_const(x, C = C) +
      f_cos(x, omega = omega_cs) +
      f_exp(x, a = a) +
      f_sin(x, omega = omega_sn)
  }
  
  
  f_C <- f_const |> generate_ts(n, C = C)
  f_c <- f_cos |> generate_ts(n, omega = omega_cs)
  f_s <- f_sin |> generate_ts(n, omega = omega_sn)
  f_e <- f_exp |> generate_ts(n, a = a)
  
  mse_lst <- list()
  for (i in 1:n) {
    f_noise <- rnorm(n, sd = sigma)
    
    f_n <- f_sum(1:n) + f_noise
    
    
    
    c <- circulant_SSA(f_n, L = L, extend_flag = TRUE)
    # r <- c$t_series
    r <- grouping_cissa(c, groups= list(trend = c(0, 1/1000), 
                                        sesonal2 = c(1/25, 1/23),
                                        sesonal1 = c(1/13, 1/10)
                                        ))$t_series
    
    # mse_lst$cissa <- c(mse_lst$cissa, mse(f_sum(1:n), r[, 9] + r[, 5] + r[, 1])) 
    mse_lst$cissa <- c(mse_lst$cissa,
                       mse(f_sum(1:n),
                           r$trend + r$sesonal1 + r$sesonal2)) 
    
    
    
    
    s <- ssa(f_n, L)
    # e <- eossa(s, 1:10, k = 6)
    e <- fossa(s)
    
    g_sesonal <- grouping.auto(e, base = "eigen",
                       freq.bins = list(trend = 1/1000, 
                                        sesonal2 = c(1/25, 1/23),
                                        sesonal1 = c(1/13, 1/10)
                                        ),
                       threshold = 0.5)
    
    r <- reconstruct(e, groups=c(list(exp = 1, C = 2), g_sesonal))
    
    mse_lst$ssa <- 
      c(mse_lst$ssa, mse(f_sum(1:n), r$trend + r$sesonal2 + r$sesonal1))
 
  }
  return(mse_lst)
}

res_mse_test <- n_mse_tests(10000)

# Оценка плотности
density_estimate_cissa <- density(res_mse_test$cissa)

# Построение графика плотности
plot(density_estimate_cissa, main = "Оценка плотности", 
     xlab = "Значение", ylab = "Плотность", 
     col = "blue", lwd = 2)

density_estimate_ssa <- density(res_mse_test$ssa)

# Построение графика плотности
plot(density_estimate_ssa, main = "Оценка плотности", 
     xlab = "Значение", ylab = "Плотность", 
     col = "blue", lwd = 2)

res_mse_test$cissa |> summary()

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.02299 0.03062 0.03345 0.03378 0.03575 0.04800 

res_mse_test$cissa |> sd()

[1] 0.004333975

res_mse_test$ssa |> summary()

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.000576 0.001710 0.002125 0.002228 0.002585 0.006311 

res_mse_test$ssa |> sd()

[1] 0.0008283841

Как выполняется расширение ряда

IP_values_slice |> extend(268) |> length()

[1] 1073

IP_values_slice |> length()

[1] 537

IP_values_slice |> extend(192) |> plot(type="l", lwd = 3)
c(rep(0, 192),IP_values_slice) |> lines(type="l", col="red")

|

n <- 96*2 + 6
x <- 0:(n-1)
y <- sin(2*pi/12 * x)
L <- 96+6
y |> extend(L) |> plot(type="l", lwd = 3)
c(rep(0, L), y) |> lines(type="l", col="red", lwd = 3)

CiSSA через Фурье

IP_values

cissa_like_fourier_transform <- function(ts, L){
  reconstruct_fft <- function(x, y, frequencies) {
    fft_y <- fft(y)
  
    amplitudes <- Mod(fft_y)
    phases <- Arg(fft_y)
    
    reconstructed <- matrix(0, length(amplitudes), length(x))
    n <- length(amplitudes)
    for (i in 1:(length(amplitudes))) {
      reconstructed[i, ] <-
        amplitudes[i] * 
        cos(2 * pi * frequencies[i] * (x) + phases[i]) /
        n
    }
    return(list(
      ts = reconstructed,
      frequencies = frequencies
      ))
  }
  
  n <- length(ts)
  x <- 0:(n-1)
  L <- L
  frequencies <- (0:(L-1)) / L
  y <- ts
  y_main <- y
  x_main <- x
  
  X <- hankel(y, L)
  K <- dim(X)[2]
  res <- list()
  for (i in 1:K){
    y <- X[, i]
    x <- x_main[1:(1 + L - 1)]
    res[[i]] <- reconstruct_fft(x, y, frequencies)$ts
  }
  
  # for (i in 1:L){
  #   plot(x, res[[i]][9, ])
  # }
  
  res_mult <- res
  # res
  
  averaging <- function(res_comp_wise_mult){
    K <- dim(X)[2]
    counters <- rep(0, n)
    res <- matrix(0, ncol = n, nrow = L)
    for (i in 1:length(res_comp_wise_mult)){
      res[, i:(i+L-1)] <- res[, i:(i+L-1)] + res_comp_wise_mult[[i]]
      counters[i:(i+L-1)] <- counters[i:(i+L-1)] + 1
    }
    for (i in 1:n){
      res[, i] <- res[, i] / counters[i]
    }
    res
  }
  
  avr <- averaging(res_mult)
  
  group_by_elementary_freq_foureir <- function(res_averaged){
    nf2 <- 0
    if (L %% 2) {
      nf2 <- (L + 1) / 2 - 1
    } else {
      nf2 <- L / 2 - 1
    }
    nft <- nf2 + abs((L %% 2) - 2)
    
    Z <- matrix(0, ncol = nft, nrow = n)
    
    # print(Z |> dim())
    # print(res_averaged |> dim())
    
    Z[, 1] <- res_averaged[1, ]
    for (k in 1:nf2) {
      Z[, k + 1] <- res_averaged[k + 1, ] + res_averaged[L + 2 - (k + 1), ]
    }
    if (L %% 2 != 0) {
      Z[, nft] <- res_averaged[nft, ]
    }
    
    
    return(list(
      t_series = Z,
      freq = (0:dim(Z)[2])/L
    ))
  }


  rs <- group_by_elementary_freq_foureir(avr)
  return(rs)
}

data_slice <- 1:537
dates_slice <- dates[data_slice]
IP_values_slice <- IP_values[data_slice]
eps <- 1/193


c <- circulant_SSA(IP_values_slice, L = 192, extend_flag = FALSE)
r <- c$t_series

c_ft <- cissa_like_fourier_transform(IP_values_slice, L = 192)
r_ft <- c_ft$t_series

for (i in 1:dim(r)[2]){
  plot(data_slice, r_ft[, i], col= "red", type = "l", lwd = 2)
  lines(data_slice, r[, i])
}

Сравнение разложений всего временного ряда через метод Фурье, SSA, auto_SSA, CiSSA и CiSSA с расширением ряда

reconstruct_fft <- function(x_init, y_init, extend_flag = FALSE) {
    x <- x_init
    y <- y_init
    N <- length(y_init)
    H <- 0
    if (extend_flag == TRUE){
      H <- length(y) %/% 2
      y <- extend(y, H)
      x <- 0:(length(y) - 1)
    }
    
    frequencies <- (0:(length(x)-1)) / length(x)
    fft_y <- fft(y)
    
    amplitudes <- Mod(fft_y)
    phases <- Arg(fft_y)
    
    reconstructed <- matrix(0, length(amplitudes), length(x))
    n <- length(amplitudes)
    L <- n
    for (i in 1:(length(amplitudes))) {
      reconstructed[i, ] <-
        amplitudes[i] * 
        cos(2 * pi * frequencies[i] * (x) + phases[i]) /
        n
    }
    
    
    
  
    group_by_elementary_freq_foureir <- function(res_averaged){
      nf2 <- 0
      if (L %% 2) {
        nf2 <- (L + 1) / 2 - 1
      } else {
        nf2 <- L / 2 - 1
      }
      nft <- nf2 + abs((L %% 2) - 2)
      
      Z <- matrix(0, ncol = nft, nrow = n)
      
      # print(Z |> dim())
      # print(res_averaged |> dim())
      
      Z[, 1] <- res_averaged[1, ]
      for (k in 1:nf2) {
        Z[, k + 1] <- res_averaged[k + 1, ] + res_averaged[L + 2 - (k + 1), ]
      }
      if (L %% 2 != 0) {
        Z[, nft] <- res_averaged[nft, ]
      }
      
      
      return(list(
        t_series = Z[(H+1):(N+H), ],
        freq = (0:(dim(Z)[2]-1))/L
      ))
    }
  
  
    rs <- group_by_elementary_freq_foureir(reconstructed)
    
    
    return(list(
      t_series = rs$t_series,
      freq = rs$freq
      ))
  }

Идеальный случай для CiSSA и Fourier

n <- 96*2
L <- 96
x <- 0:(n-1)
y1 <- sin(2*pi/12 * x)
y2 <- cos(2*pi/3 * x)/2
y <- y1 + y2
X <- hankel(y, L = L)
eps <- 1/(n+1)

s_ssa <- ssa(y[1:(n-1)], L)
r_ssa <- reconstruct(s_ssa, groups = list(
  F1 = c(1, 2),
  F2 = c(3, 4)
))
# r_ssa <- reconstruct(s_ssa, groups=list(
#   sesonal_sin = c(1, 2),
#   sesonal_cos = c(3, 4)
# ))
# plot(x, r_ssa$F1, type="l")
# plot(x, r_ssa$F2, type= "l")
# plot(wcor(s_ssa, groups = 1:10), scales = list(at = c(10, 20, 30)))
e_ssa <- eossa_new(s_ssa, nested.groups = list(1:4), clust_type = "distance")
g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                       freq.bins = list(
                         sesonal_sin = c(1/12-eps, 1/12+eps),
                         sesonal_cos = c(1/3-eps, 1/3+eps)
                                        ),
                       threshold = 0.5)
r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
# plot(x, r_ssa_e$sesonal_sin)
# plot(x, r_ssa_e$sesonal_cos)



r_fft <- reconstruct_fft(x, y)
r_fft_grouped <- grouping_cissa(r_fft,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

r_fft_extended <- reconstruct_fft(x, y, TRUE)
r_fft_grouped_extended <- grouping_cissa(r_fft_extended,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# r_fft_grouped$sesonal_sin |> length()

# plot(x, r_fft_grouped$sesonal_sin)
# plot(x, r_fft_grouped$sesonal_cos)
# plot(x, r_fft_grouped$residuals)


r_cissa <- circulant_SSA(y, L)
r_cissa_grouped <- grouping_cissa(r_cissa,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series


r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# plot(x, r_cissa_grouped$sesonal_sin)
# plot(x, r_cissa_grouped$sesonal_cos)
# plot(x, r_cissa_grouped$residuals)

library(xtable)

# Шаг 2: Создание примера данных
data <- data.frame(
  Метод = c("SSA", "SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  sin_err = c (1, 20, 20, 1, 1, 1),
  cos_err = c(1, 1, 1, 1, 1, 1)
)

data$sin_err[1] <- mse(y1[1:(n-1)], r_ssa$F1) |>
  formatC(format = "e", digits = 1) 
data$cos_err[1] <- mse(y2[1:(n-1)], r_ssa$F2) |>
  formatC(format = "e", digits = 1)

data$cos_err[2] <- mse(y1[1:(n-1)], r_ssa_e$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mse(y2[1:(n-1)], r_ssa_e$sesonal_cos) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mse(y1, r_fft_grouped$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mse(y2, r_fft_grouped$sesonal_cos) |>
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mse(y1, r_cissa_grouped$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mse(y2, r_cissa_grouped$sesonal_cos) |> 
  formatC(format = "e", digits = 1)


data$cos_err[5] <- mse(y1, r_cissa_grouped_ext$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[5] <- mse(y2, r_cissa_grouped_ext$sesonal_cos) |> 
  formatC(format = "e", digits = 1)

data$cos_err[6] <- mse(y1, r_fft_grouped_extended$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[6] <- mse(y2, r_fft_grouped_extended$sesonal_cos) |>
  formatC(format = "e", digits = 1)



table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Thu Nov 28 19:49:16 2024
\begin{table}[ht]
\centering
\begin{tabular}{lll}
  \hline
Метод & sin\_err & cos\_err \\ 
  \hline
SSA & 6.8e-30 & 1.5e-29 \\ 
  SSA EOSSA & 1.5e-29 & 7.5e-30 \\ 
  Fourier & 1.7e-28 & 3.5e-28 \\ 
  CiSSA & 1.9e-29 & 5.3e-30 \\ 
  CiSSA extended & 2.0e-04 & 8.6e-04 \\ 
  Fourier extended & 6.2e-04 & 2.6e-03 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Идеальный случай с шумом

data <- list(
  Метод = c("SSA","SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  sin_err = list(list(), list(), list(), list(), list(), list()),
  cos_err = list(list(), list(), list(), list(), list(), list())
)

for (i in 1:50){
  set.seed(i)
  
  n <- 96*2
  L <- 96
  x <- 0:(n-1)
  y1 <- sin(2*pi/12 * x)
  y2 <- cos(2*pi/3 * x)/2
  y <- y1 + y2 + rnorm(n, 0, 0.1)
  X <- hankel(y, L = L)
  eps <- 1/(n+1)
  
  s_ssa <- ssa(y[1:(n-1)], L)
  r_ssa <- reconstruct(s_ssa, groups = list(
    F1 = c(1, 2),
    F2 = c(3, 4)
  ))
  # r_ssa <- reconstruct(s_ssa, groups=list(
  #   sesonal_sin = c(1, 2),
  #   sesonal_cos = c(3, 4)
  # ))
  # plot(x[1:(n-1)], r_ssa$F1, type = "l")
  # plot(x[1:(n-1)], r_ssa$F2, type = "l")
  # plot(wcor(s_ssa, groups = 1:10), scales = list(at = c(10, 20, 30)))
  e_ssa <- eossa_new(s_ssa, nested.groups = list(1:4), clust_type = "distance")
  g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                         freq.bins = list(
                           sesonal_sin = c(1/12-eps, 1/12+eps),
                           sesonal_cos = c(1/3-eps, 1/3+eps)
                                          ),
                         threshold = 0.5)
  r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
  # plot(x, r_ssa_e$sesonal_sin)
  # plot(x, r_ssa_e$sesonal_cos)
  
  
  
  r_fft <- reconstruct_fft(x, y)
  r_fft_grouped <- grouping_cissa(r_fft,
                      groups = list(
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  r_fft_extended <- reconstruct_fft(x, y, TRUE)
  r_fft_grouped_extended <- grouping_cissa(r_fft_extended,
                      groups = list(
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  # r_fft_grouped$sesonal_sin |> length()
  
  # plot(x, r_fft_grouped$sesonal_sin)
  # plot(x, r_fft_grouped$sesonal_cos)
  # plot(x, r_fft_grouped$residuals)
  
  
  r_cissa <- circulant_SSA(y, L)
  r_cissa_grouped <- grouping_cissa(r_cissa,
                      groups = list(
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  
  r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
  r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                      groups = list(
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  data$cos_err[[1]][[i]] <- 
    min(mse(y1[1:(n-1)], r_ssa$F1), mse(y1[1:(n-1)], r_ssa$F2))
  data$sin_err[[1]][[i]] <-
    min(mse(y2[1:(n-1)], r_ssa$F1), mse(y2[1:(n-1)], r_ssa$F2))
  
  data$cos_err[[2]][[i]] <- mse(y1[1:(n-1)], r_ssa_e$sesonal_sin)
  data$sin_err[[2]][[i]] <- mse(y2[1:(n-1)], r_ssa_e$sesonal_cos)
  
  
  data$cos_err[[3]][[i]] <- mse(y1, r_fft_grouped$sesonal_sin)
  data$sin_err[[3]][[i]] <- mse(y2, r_fft_grouped$sesonal_cos)
  
  
  data$cos_err[[4]][[i]] <- mse(y1, r_cissa_grouped$sesonal_sin)
  data$sin_err[[4]][[i]] <- mse(y2, r_cissa_grouped$sesonal_cos)
  
  
  data$cos_err[[5]][[i]] <- mse(y1, r_cissa_grouped_ext$sesonal_sin)
  data$sin_err[[5]][[i]] <- mse(y2, r_cissa_grouped_ext$sesonal_cos)
  
  data$cos_err[[6]][[i]] <- mse(y1, r_fft_grouped_extended$sesonal_sin)
  data$sin_err[[6]][[i]] <- mse(y2, r_fft_grouped_extended$sesonal_cos)
}
library(xtable)

# Шаг 2: Создание примера данных
data_prev <- data

data <- data.frame(
  Метод = c("SSA", "SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  sin_err = c(0, 0, 0, 0, 0, 0),
  cos_err = c(0, 0, 0, 0, 0, 0)
)

data$cos_err[1] <- mean(data_prev$cos_err[[1]] |> unlist()) |>
  formatC(format = "e", digits = 1)
data$sin_err[1] <- mean(data_prev$sin_err[[1]]|> unlist()) |>
  formatC(format = "e", digits = 1)

data$cos_err[2] <- mean(data_prev$cos_err[[2]] |> unlist()) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mean(data_prev$sin_err[[2]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mean(data_prev$cos_err[[3]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mean(data_prev$sin_err[[3]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mean(data_prev$cos_err[[4]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mean(data_prev$sin_err[[4]]|> unlist()) |> 
  formatC(format = "e", digits = 1)


data$cos_err[5] <- mean(data_prev$cos_err[[5]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[5] <- mean(data_prev$sin_err[[5]]|> unlist()) |> 
  formatC(format = "e", digits = 1)

data$cos_err[6] <- mean(data_prev$cos_err[[6]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[6] <- mean(data_prev$sin_err[[6]]|> unlist()) |> 
  formatC(format = "e", digits = 1)

for (i in 1:4){
  for (j in (i):4){
    x <- data_prev$cos_err[[i]] |> unlist()
    y <- data_prev$cos_err[[j]] |> unlist()
    t_test_result <- t.test(x, y, paired = TRUE)
    print(paste("cos, ", data$Метод[i], " ", data$Метод[j], ", p-val = ", t_test_result$p.value))
  }
}

[1] "cos,  SSA   SSA , p-val =  NaN"
[1] "cos,  SSA   SSA EOSSA , p-val =  0.038881875174499"
[1] "cos,  SSA   Fourier , p-val =  3.63996112806195e-10"
[1] "cos,  SSA   CiSSA , p-val =  1.69300523929475e-11"
[1] "cos,  SSA EOSSA   SSA EOSSA , p-val =  NaN"
[1] "cos,  SSA EOSSA   Fourier , p-val =  3.61837650659192e-10"
[1] "cos,  SSA EOSSA   CiSSA , p-val =  1.68793475607746e-11"
[1] "cos,  Fourier   Fourier , p-val =  NaN"
[1] "cos,  Fourier   CiSSA , p-val =  1.84590409706372e-05"
[1] "cos,  CiSSA   CiSSA , p-val =  NaN"

for (i in 1:4){
  for (j in (i):4){
    x <- data_prev$sin_err[[i]] |> unlist()
    y <- data_prev$sin_err[[j]] |> unlist()
    t_test_result <- t.test(x, y, paired = TRUE)
    print(paste("sin, ", data$Метод[i], " ", data$Метод[j], ", p-val = ", t_test_result$p.value))
  }
}

[1] "sin,  SSA   SSA , p-val =  NaN"
[1] "sin,  SSA   SSA EOSSA , p-val =  3.71816430587151e-11"
[1] "sin,  SSA   Fourier , p-val =  5.97385946308374e-13"
[1] "sin,  SSA   CiSSA , p-val =  2.11234626356184e-13"
[1] "sin,  SSA EOSSA   SSA EOSSA , p-val =  NaN"
[1] "sin,  SSA EOSSA   Fourier , p-val =  5.89508148968624e-13"
[1] "sin,  SSA EOSSA   CiSSA , p-val =  2.13061715482733e-13"
[1] "sin,  Fourier   Fourier , p-val =  NaN"
[1] "sin,  Fourier   CiSSA , p-val =  2.79473942754082e-06"
[1] "sin,  CiSSA   CiSSA , p-val =  NaN"

table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Thu Nov 28 19:53:46 2024
\begin{table}[ht]
\centering
\begin{tabular}{lll}
  \hline
Метод & sin\_err & cos\_err \\ 
  \hline
SSA & 3.2e-04 & 3.0e-04 \\ 
  SSA EOSSA & 3.2e-04 & 3.0e-04 \\ 
  Fourier & 1.3e-04 & 1.1e-04 \\ 
  CiSSA & 1.8e-04 & 1.8e-04 \\ 
  CiSSA extended & 7.1e-04 & 1.9e-03 \\ 
  Fourier extended & 1.3e-03 & 3.9e-03 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Добавим непериодичность

n <- 96*2
x <- 0:(n-1)
L <- 96
y1 <- sin(2*pi/12 * x)
y2 <- cos(2*pi/3 * x)/2
y3 <- exp(x/100) + 1
y <- y1 + y2 + y3
eps <- 1/(n+1)

s_ssa <- ssa(y[1:(n-1)], L)
r_ssa <- reconstruct(s_ssa, groups=list(
  e = c(1),
  sesonal_sin = c(2, 3),
  sesonal_cos = c(4, 5)
))
# plot(x[1:(n-1)], r_ssa$e, type = "l")
# plot(x[1:(n-1)], r_ssa$sesonal_sin, type = "l")
# plot(x[1:(n-1)], r_ssa$sesonal_cos, type = "l")
# plot(wcor(s_ssa, groups = 1:10), scales = list(at = c(10, 20, 30)))
e_ssa <- eossa_new(s_ssa, nested.groups = list(1:7), clust_type = "distance")
g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                       freq.bins = list(
                         e = c(0, 1/12-eps-eps),
                         sesonal_sin = c(1/12-eps, 1/12+eps),
                         sesonal_cos = c(1/3-eps, 1/3+eps)
                                        ),
                       threshold = 0.5)
r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
# plot(x, r_ssa_e$sesonal_sin)
# plot(x, r_ssa_e$sesonal_cos)



r_fft <- reconstruct_fft(x, y)
r_fft_grouped <- grouping_cissa(r_fft,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

r_fft_extended <- reconstruct_fft(x, y, TRUE)
r_fft_grouped_extended <- grouping_cissa(r_fft_extended,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# r_fft_grouped$sesonal_sin |> length()

# plot(x, r_fft_grouped$sesonal_sin)
# plot(x, r_fft_grouped$sesonal_cos)
# plot(x, r_fft_grouped$residuals)


r_cissa <- circulant_SSA(y, L)
r_cissa_grouped <- grouping_cissa(r_cissa,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series


r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# plot(x, r_cissa_grouped$sesonal_sin)
# plot(x, r_cissa_grouped$sesonal_cos)
# plot(x, r_cissa_grouped$residuals)

library(xtable)

# Шаг 2: Создание примера данных
data <- data.frame(
  Метод = c("SSA", "SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  e_err = c(1, 1, 1, 1, 1, 1),
  sin_err = c (1, 20, 20, 1, 1, 1),
  cos_err = c(1, 1, 1, 1, 1, 1)
)

data$cos_err[1] <- mse(y1[1:(n-1)], r_ssa$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[1] <- mse(y2[1:(n-1)], r_ssa$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$e_err[1] <- mse(y3[1:(n-1)], r_ssa$e) |>
  formatC(format = "e", digits = 1)

data$cos_err[2] <- mse(y1[1:(n-1)], r_ssa_e$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mse(y2[1:(n-1)], r_ssa_e$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$e_err[2] <- mse(y3[1:(n-1)], r_ssa_e$e) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mse(y1, r_fft_grouped$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mse(y2, r_fft_grouped$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$e_err[3] <- mse(y3, r_fft_grouped$e) |>
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mse(y1, r_cissa_grouped$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mse(y2, r_cissa_grouped$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$e_err[4] <- mse(y3, r_cissa_grouped$e) |>
  formatC(format = "e", digits = 1)

data$cos_err[5] <- mse(y1, r_cissa_grouped_ext$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[5] <- mse(y2, r_cissa_grouped_ext$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$e_err[5] <- mse(y3, r_cissa_grouped_ext$e) |>
  formatC(format = "e", digits = 1)

data$cos_err[6] <- mse(y1, r_fft_grouped_extended$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[6] <- mse(y2, r_fft_grouped_extended$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$e_err[6] <- mse(y3, r_fft_grouped_extended$e) |>
  formatC(format = "e", digits = 1)



table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Thu Nov 28 19:55:50 2024
\begin{table}[ht]
\centering
\begin{tabular}{llll}
  \hline
Метод & e\_err & sin\_err & cos\_err \\ 
  \hline
SSA & 5.0e-03 & 8.9e-07 & 5.2e-05 \\ 
  SSA EOSSA & 1.7e-28 & 1.6e-29 & 8.7e-30 \\ 
  Fourier & 1.1e-01 & 6.1e-04 & 6.8e-03 \\ 
  CiSSA & 5.3e-02 & 1.6e-05 & 4.9e-04 \\ 
  CiSSA extended & 5.0e-04 & 2.1e-04 & 1.1e-03 \\ 
  Fourier extended & 1.4e-03 & 1.3e-03 & 8.4e-03 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Непериодичность с шумом

# data <- data.frame(
#   Метод = c("Fourier", "CiSSA", "CiSSA с расширением ряда"),
#   sin_err = c (20, 20, 20),
#   cos_err = c(1, 1, 20),
#   exp_err = c(1, 1, 20)
# )
data <- list(
  Метод = c("SSA", "SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  sin_err = list(list(), list(), list(), list(), list(), list()),
  cos_err = list(list(), list(), list(), list(), list(), list()),
  exp_err = list(list(), list(), list(), list(), list(), list())
)

for (i in 1:100){
  set.seed(i)
  
  n <- 96*2
  x <- 0:(n-1)
  L <- 96
  y1 <- sin(2*pi/12 * x)
  y2 <- cos(2*pi/3 * x)/2
  y3 <- exp(x/100) + 1
  y <- y1 + y2 + y3 + rnorm(n, 0, 0.1)
  eps <- 1/(n+1)
  
  s_ssa <- ssa(y[1:(n-1)], L)
  r_ssa <- reconstruct(s_ssa, groups=list(
    e = 1,
    sesonal_sin = c(2, 3),
    sesonal_cos = c(4, 5)
  ))
  
  
  e_ssa <- eossa_new(s_ssa, nested.groups = list(1:6), clust_type = "distance")
  g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                         freq.bins = list(
                           e = c(0, 1/12-eps-eps),
                           sesonal_sin = c(1/12-eps, 1/12+eps),
                           sesonal_cos = c(1/3-eps, 1/3+eps)
                                          ),
                         threshold = 0.5)
  r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
  
  
  
  r_fft <- reconstruct_fft(x, y)
  r_fft_grouped <- grouping_cissa(r_fft,
                      groups = list(
                        e = c(0, 1/12-eps-eps),
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  r_fft_extended <- reconstruct_fft(x, y, TRUE)
  r_fft_grouped_extended <- grouping_cissa(r_fft_extended,
                      groups = list(
                        e = c(0, 1/12-eps-eps),
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  
  r_cissa <- circulant_SSA(y, L)
  r_cissa_grouped <- grouping_cissa(r_cissa,
                      groups = list(
                        e = c(0, 1/12-eps-eps),
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  
  r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
  r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                      groups = list(
                        e = c(0, 1/12-eps-eps),
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  
  
  data$cos_err[[1]][[i]] <- mse(y1[1:(n-1)], r_ssa$sesonal_sin)
  data$sin_err[[1]][[i]] <- mse(y2[1:(n-1)], r_ssa$sesonal_cos)
  data$exp_err[[1]][[i]] <- mse(y3[1:(n-1)], r_ssa$e)
  
  data$cos_err[[2]][[i]] <- mse(y1[1:(n-1)], r_ssa_e$sesonal_sin)
  data$sin_err[[2]][[i]] <- mse(y2[1:(n-1)], r_ssa_e$sesonal_cos)
  data$exp_err[[2]][[i]] <- mse(y3[1:(n-1)], r_ssa_e$e)
  
  # print(data$sin_err[[2]][[i]])
  
  
  data$cos_err[[3]][[i]] <- mse(y1, r_fft_grouped$sesonal_sin)
  data$sin_err[[3]][[i]] <- mse(y2, r_fft_grouped$sesonal_cos)
  data$exp_err[[3]][[i]] <- mse(y3, r_fft_grouped$e)
  
  
  data$cos_err[[4]][[i]] <- mse(y1, r_cissa_grouped$sesonal_sin)
  data$sin_err[[4]][[i]] <- mse(y2, r_cissa_grouped$sesonal_cos)
  data$exp_err[[4]][[i]] <- mse(y3, r_cissa_grouped$e)
  
  
  data$cos_err[[5]][[i]] <- mse(y1, r_cissa_grouped_ext$sesonal_sin)
  data$sin_err[[5]][[i]] <- mse(y2, r_cissa_grouped_ext$sesonal_cos)
  data$exp_err[[5]][[i]] <- mse(y3, r_cissa_grouped_ext$e)
  
  data$cos_err[[6]][[i]] <- mse(y1, r_fft_grouped_extended$sesonal_sin)
  data$sin_err[[6]][[i]] <- mse(y2, r_fft_grouped_extended$sesonal_cos)
  data$exp_err[[6]][[i]] <- mse(y3, r_fft_grouped_extended$e)
}
library(xtable)
# data$sin_err[[2]]
# Шаг 2: Создание примера данных
data_prev <- data
data <- data.frame(
  Метод = c("SSA", "SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  sin_err = c(0, 0, 0, 0, 0, 0),
  cos_err = c(0, 0, 0, 0, 0, 0),
  exp_err = c(0, 0, 0, 0, 0, 0)
)

data$cos_err[1] <- mean(data_prev$cos_err[[1]] |> unlist()) |>
  formatC(format = "e", digits = 1)
data$sin_err[1] <- mean(data_prev$sin_err[[1]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$exp_err[1] <- mean(data_prev$exp_err[[1]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[2] <- mean(data_prev$cos_err[[2]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mean(data_prev$sin_err[[2]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$exp_err[2] <- mean(data_prev$exp_err[[2]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mean(data_prev$cos_err[[3]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mean(data_prev$sin_err[[3]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$exp_err[3] <- mean(data_prev$exp_err[[3]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mean(data_prev$cos_err[[4]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mean(data_prev$sin_err[[4]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$exp_err[4] <- mean(data_prev$exp_err[[4]]|> unlist()) |>
  formatC(format = "e", digits = 1)

data$cos_err[5] <- mean(data_prev$cos_err[[5]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[5] <- mean(data_prev$sin_err[[5]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$exp_err[5] <- mean(data_prev$exp_err[[5]]|> unlist()) |>
  formatC(format = "e", digits = 1)

data$cos_err[6] <- mean(data_prev$cos_err[[6]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[6] <- mean(data_prev$sin_err[[6]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$exp_err[6] <- mean(data_prev$exp_err[[6]]|> unlist()) |>
  formatC(format = "e", digits = 1)

for (i in 1:4){
  for (j in (i):4){
    x <- data_prev$cos_err[[i]] |> unlist()
    y <- data_prev$cos_err[[j]] |> unlist()
    t_test_result <- t.test(x, y, paired = TRUE)
    print(paste("cos, ", data$Метод[i], " ", data$Метод[j], ", p-val = ", t_test_result$p.value))
  }
}

[1] "cos,  SSA   SSA , p-val =  NaN"
[1] "cos,  SSA   SSA EOSSA , p-val =  1.15068669354673e-13"
[1] "cos,  SSA   Fourier , p-val =  1.15081927627114e-78"
[1] "cos,  SSA   CiSSA , p-val =  4.65220864068305e-44"
[1] "cos,  SSA EOSSA   SSA EOSSA , p-val =  NaN"
[1] "cos,  SSA EOSSA   Fourier , p-val =  6.8490552789737e-79"
[1] "cos,  SSA EOSSA   CiSSA , p-val =  7.58217507825911e-49"
[1] "cos,  Fourier   Fourier , p-val =  NaN"
[1] "cos,  Fourier   CiSSA , p-val =  4.1470816136733e-79"
[1] "cos,  CiSSA   CiSSA , p-val =  NaN"

for (i in 1:4){
  for (j in (i):4){
    x <- data_prev$sin_err[[i]] |> unlist()
    y <- data_prev$sin_err[[j]] |> unlist()
    t_test_result <- t.test(x, y, paired = TRUE)
    print(paste("sin, ", data$Метод[i], " ", data$Метод[j], ", p-val = ", t_test_result$p.value))
  }
}

[1] "sin,  SSA   SSA , p-val =  NaN"
[1] "sin,  SSA   SSA EOSSA , p-val =  0.273046993097837"
[1] "sin,  SSA   Fourier , p-val =  3.41671346310886e-19"
[1] "sin,  SSA   CiSSA , p-val =  2.68804866311398e-16"
[1] "sin,  SSA EOSSA   SSA EOSSA , p-val =  NaN"
[1] "sin,  SSA EOSSA   Fourier , p-val =  3.14261131207794e-19"
[1] "sin,  SSA EOSSA   CiSSA , p-val =  5.94062534986952e-16"
[1] "sin,  Fourier   Fourier , p-val =  NaN"
[1] "sin,  Fourier   CiSSA , p-val =  2.93671115450334e-28"
[1] "sin,  CiSSA   CiSSA , p-val =  NaN"

for (i in 1:4){
  for (j in (i):4){
    x <- data_prev$exp_err[[i]] |> unlist()
    y <- data_prev$exp_err[[j]] |> unlist()
    t_test_result <- t.test(x, y, paired = TRUE)
    print(paste("exp, ", data$Метод[i], " ", data$Метод[j], ", p-val = ", t_test_result$p.value))
  }
}

[1] "exp,  SSA   SSA , p-val =  NaN"
[1] "exp,  SSA   SSA EOSSA , p-val =  1.2968052853879e-78"
[1] "exp,  SSA   Fourier , p-val =  3.1975845882856e-218"
[1] "exp,  SSA   CiSSA , p-val =  6.62421968424602e-130"
[1] "exp,  SSA EOSSA   SSA EOSSA , p-val =  NaN"
[1] "exp,  SSA EOSSA   Fourier , p-val =  5.2539864487142e-262"
[1] "exp,  SSA EOSSA   CiSSA , p-val =  2.10736945563773e-133"
[1] "exp,  Fourier   Fourier , p-val =  NaN"
[1] "exp,  Fourier   CiSSA , p-val =  4.58655868843589e-140"
[1] "exp,  CiSSA   CiSSA , p-val =  NaN"

table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Thu Nov 28 20:01:06 2024
\begin{table}[ht]
\centering
\begin{tabular}{llll}
  \hline
Метод & sin\_err & cos\_err & exp\_err \\ 
  \hline
SSA & 2.9e-04 & 3.6e-04 & 5.2e-03 \\ 
  SSA EOSSA & 2.9e-04 & 3.1e-04 & 9.4e-04 \\ 
  Fourier & 6.9e-04 & 7.2e-03 & 1.2e-01 \\ 
  CiSSA & 1.7e-04 & 7.0e-04 & 5.5e-02 \\ 
  CiSSA extended & 6.8e-04 & 2.1e-03 & 2.7e-03 \\ 
  Fourier extended & 1.9e-03 & 9.6e-03 & 3.0e-03 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Пример, который не подходит ни к какому методу

n <- 96*2 + 7
L <- 89
x <- 0:(n-1)
y1 <- sin(2*pi/13 * x)
y2 <- cos(2*pi/8 * x)
y3 <- -x*x/1000
y4 <- exp(x/55)
y <- y1 + y2 + y3 + y4 
plot(x, y)

X <- hankel(y, L = L)
eps <- 1/(n+1)

s_ssa <- ssa(y, L)
# r_ssa <- reconstruct(s_ssa, groups=list(
#   sesonal_sin = c(1, 2),
#   sesonal_cos = c(3, 4)
# ))
# plot(x, r_ssa$sesonal_sin)
# plot(x, r_ssa$sesonal_cos)
# plot(wcor(s_ssa, groups = 1:10), scales = list(at = c(10, 20, 30)))
e_ssa <- eossa_new(s_ssa, nested.groups = list(1:10), clust_type = "distance")
g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                       freq.bins = list(
                         sesonal_sin = c(1/12-eps, 1/12+eps),
                         sesonal_cos = c(1/8-eps, 1/8+eps)
                                        ),
                       threshold = 0.5)
r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
# plot(x, r_ssa_e$sesonal_sin)
# plot(x, r_ssa_e$sesonal_cos)



r_fft <- reconstruct_fft(x, y)
r_fft_grouped <- grouping_cissa(r_fft,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/8-eps, 1/8+eps)
                    )
                    )$t_series

# r_fft_grouped$sesonal_sin |> length()

# plot(x, r_fft_grouped$sesonal_sin)
# plot(x, r_fft_grouped$sesonal_cos)
# plot(x, r_fft_grouped$residuals)


r_cissa <- circulant_SSA(y, L)
r_cissa_grouped <- grouping_cissa(r_cissa,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/8-eps, 1/8+eps)
                    )
                    )$t_series


r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/8-eps, 1/8+eps)
                    )
                    )$t_series

# plot(x, r_cissa_grouped$sesonal_sin)
# plot(x, r_cissa_grouped$sesonal_cos)
# plot(x, r_cissa_grouped$residuals)

library(xtable)

# Шаг 2: Создание примера данных
data <- data.frame(
  Метод = c("SSA EOSSA","Fourier", "CiSSA", "CiSSA extended"),
  sin_err = c (20, 20, 1, 1),
  cos_err = c(1, 1, 1, 1)
)

data$cos_err[1] <- mse(y1, r_ssa_e$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[1] <- mse(y2, r_ssa_e$sesonal_cos) |>
  formatC(format = "e", digits = 1)


data$cos_err[2] <- mse(y1, r_fft_grouped$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mse(y2, r_fft_grouped$sesonal_cos) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mse(y1, r_cissa_grouped$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mse(y2, r_cissa_grouped$sesonal_cos) |> 
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mse(y1, r_cissa_grouped_ext$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mse(y2, r_cissa_grouped_ext$sesonal_cos) |> 
  formatC(format = "e", digits = 1)



table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Thu Nov 28 18:27:55 2024
\begin{table}[ht]
\centering
\begin{tabular}{lll}
  \hline
Метод & sin\_err & cos\_err \\ 
  \hline
SSA EOSSA & 2.1e-21 & 7.4e-22 \\ 
  Fourier & 1.6e-02 & 4.2e-01 \\ 
  CiSSA & 2.2e-02 & 3.7e-02 \\ 
  CiSSA extended & 2.9e-03 & 4.5e-03 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}